Role of Galactic sources and magnetic fields in forming the observed 
energy-dependent composition of ultrahigh-energy cosmic rays 
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Recent results from the Pierre Auger Observatory, showing energy-dependent chemical composi- 
tion of ultrahigh-energy cosmic rays (UHECRs) with a growing fraction of heavy elements at high 
energies, suggest a possible non-negligible contribution of the Galactic sources. We show that, in 
the case of UHECRs produced by gamma-ray bursts (GRBs) or rare types of supernova explosions 
that took place in the Milky Way in the past, the change in UHECR composition can result from 
the difference in diffusion times for different species. The anisotropy in the direction of the Galactic 
Center is expected to be a few per cent on average, but the locations of the most recent/closest 
bursts can be associated with observed clusters of UHECRs. 
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Recent measurements of the composition of ultrahigh- 
energy cosmic rays (UHECRs) by the Pierre Auger Ob- 
servatory (PAO) have suggested that the mean nuclear 
mass may increase with energy between 2 EeV and 
35 EeV [1]. This was not expected. The lack of plausible 
sources in the Milky Way (MW) and the lack of Galac- 
tocentric anisotropy of the arrival directions of UHECRs 
are usually seen as evidence for extragalactic origin of 
UHECRs above 10 18 eV. 

However, if the cosmic rays of 10 18 — 10 19 eV are nuclei, 
the turbulent Galactic micro-Gauss magnetic fields @ 
can retain them in the Galaxy and isotropize their di- 
rections sufficiently to show no disagreement with the 
data. Moreover, since the diffusion times depend on the 
rigidity, the observed composition can be altered by dif- 
fusion Q. Since the heavier nuclei spend more time in 
the Galaxy than the lighter nuclei and protons, they have 
the higher number density and flux. Thus, diffusion alone 
can alter the composition of UHECRs produced by the 
Galactic sources and increase the observed fraction of 
nuclei. 

As for the plausible sources, there is a growing evi- 
dence that long GRBs are caused by a relatively rare 
type of supernovae, while the short GRBs probably re- 
sult from the coalescence of neutron stars and/or black 
holes [4]. Compact star mergers undoubtedly take place 
in the Milky Way, and therefore short GRBs should oc- 
cur in our Galaxy. Although there is some correlation 
of long GRBs with star- forming metal-poor galaxies Q , 
many long GRBs are observed in high-mctallicity galax- 
ies as well |fj| , and therefore one expects that long GRBs 
should occur in the Milky Way. Less powerful hyper- 
novae, too weak to produce a GRB, but can still acceler- 
ate UHECRs [7| , with a substantial fraction of nuclei [8| . 

We will show that the observed change in chemical 
composition can be explained, at least qualitatively, by 
the past Galactic GRBs producing UHECRs that diffuse 
in the ~ 3/iG magnetic field. Changes in composition 
due to a magnetic fields have been discussed in connec- 



tion with the spectral "knee" [3j , and also for a transient 
source We illustrate the effect of Galactic diffusion 
using a simple model and show that, in some range of 
energies, the observed composition is energy-dependent 
even if the spectrum produced at the source does not have 
an energy-dependent composition. Moreover, for reason- 
able values of the Galactic magnetic fields, the transi- 
tion to energy-dependent composition is expected around 
~ 10 18 eV, in agreement with the data. We concen- 
trate on UHECRs with the energies below the Greisen- 
Zatsepin-Kuzmin (GZK) cutoff jlfj ]. We note that the 
PAO results on chemical composition [l[ are in agree- 
ment with the Yakutsk experiment fllj , but they have 
not been corroborated by the HiRes experiment [12[ , and 
there remains a significant experimental uncertainty due 
to possible systematic errors. 

GRBs have been proposed as the sources of extragalac- 
tic UHECRs [l3j], and theyhave also been considered as 
possible Galactic sources [14]]. It is believed that GRBs, 
hypernovae, or other stellar events capable of produc- 
ing UHECRs could have happened in the Milky Way at 
the rate of one per tcRB ~ 10 4 — 10 5 years [15]. Such 
events have been linked to the observations of INTE- 
GRAL, Fermi and PAMELA 

Let us consider a simple model of cosmic ray transport. 
The spectrum of UHECRs is assumed to contain differ- 
ent nuclei, which we will label by their electric charges 
qi = eZi. Let us denote the number density of i'th species 
of cosmic rays with energy E at point r by rii(E, r, t). We 
assume for simplicity that production of cosmic rays of 
different species occurs in such a way that the produc- 
tion spectra Qi(E,r) have the same energy dependence, 
and that the chemical composition is energy independent, 
given by the constants £j : 

Q i (E,?) = Q Z i (E/E )-ip(f f ), (1) 

where Qo and Eg do not depend on the type of the species 
or the energy. We will consider three different models for 
the distribution of sources pif). 
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We now ask whether the measured fluxes inside the 
galaxy can be altered by diffusion. In diffusive approxi- 
mation, the transport inside the galaxy can be described 
by the usual equation: 



V(AVni) + — (km) = 
Qi(E,f,t) + J2 J P ik (E,E')n k (E')dE'. 



(2) 



Here Di(E,r,t) = Di(E) is the diffusion coefficient, 
which we will assume to be constant in space and time. 
The energy losses and all the interactions that change 
the particle energies are given by bi (E) and the kernel in 
the collision integral Pik(E, E 1 ). For energies below GZK 
cutoff, one can neglect the energy losses on the diffusion 
time scales. 

The diffusion coefficient D(E) depends primarily on 
the structure of the magnetic fields in the galaxy. Let 
us assume that the magnetic field structure is comprised 
of uniform randomly oriented domains of radius Iq with 
a constant field B in each domain. The density of such 
domains is N < l^ 3 . The Larmor radius depends on the 
particle energy E and its electric charge qi = eZi\ 



Ri 



E 



la 



E 
Eqa 



E () = Bin = 10 18 eV 



, where Eq 
B 



3 x I0- 6 GJ V 0.3 kpc 



EoZi, (3) 



(4) 



The spatial energy spectrum of random magnetic fields 
inferred from observations suggests that B ~ 3/xG on 
the 0.3 kpc spatial scales, and that there is a significant 
change at I = I/fc ~ 0.1 — 0.5 kpc 0. This can be 
understood theoretically because the turbulent energy is 
injected into the interstellar medium by supernova explo- 
sions on the scales of order 0.1 kpc. This energy is trans- 
ferred to smaller scales by direct cascade, and to larger 
scales by inverse cascade of magnetic helicity. Single-cell- 
size models favor ~ 0.1 kpc scales as well 

For each species, there is a critical energy Eo,i for which 
the Larmor radius is equal to the magnetic coherence 
length In. In this simple model, the diffusion proceeds in 
two very different regimes. For E < Eo,%, the mean free 
path of the diffusing particle is I ~ l/(iVZ§) ~ In, and, in 
the units for which the speed of light c = 1, the diffusion 
coefficient is 



Di(E) 



1 



-In 



Dn, for E <C E 0ti 



(5) 



For E ^> En t i, the particle is deflected only by a small 
angle 9 ~ lo/Ri, and, after k deflections, the mean deflec- 
tion angle squared is 9 2 ~ k(ln/ ' Ri) 2 . The corresponding 
diffusion coefficient is 



Di(E) = D 



E 

Eqa 



, for E > E . 



(6) 
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FIG. 1: Predicted UHECR spectra, assuming that the sources 
produce 90% protons and 10% iron, with identicai spectra 
oc E~ 2 ' 3 , and that the source distribution traces the distribu- 
tion of stars in the Gaiaxy. We used samples of 10 3 GRBs at 
random locations with time intervais of 10 years. The mag- 
netic fiefd was assumed to be 4 fiG, coherent over In = 0.2 kpc 
domains, 5i = 0.3, 62 = 0. The overaff power and the iron 
fraction were adjusted to fit the PAO data points [HJ (shown) . 
For each random sampfe, the fit parameters differ slightly, de- 
pending on the location of the latest/closest burst. To model 
the spectrum at E > 3 x 10 19 eV, one has to account for en- 
ergy losses and (proton) contribution of extragalactic sources, 
which we leave for future work. 



While this simple model is reasonable and well 
known 17], in reality the energy dependence of the dif- 
fusion coefficient can be more complicated, for exam- 
ple, because the magnetic field domains are different in 
size @. To account for this, we introduce two parameters 
< 81.2 < 0.5 and assume the following behavior of the 
diffusion coefficient: 



Di(E) = 



D, 



D, 



(*) • 
»(*) 



Si 

(2-s 2 ) 



E < En^i, 

E > Eqa- 



(7) 



It is essential for our discussion that the energy depen- 
dence of the diffusion coefficient changes at a critical en- 
ergy Eqa = eZiE which is different for different nuclei. 

For ~ 0.3 kpc, the mean time of diffusion of nuclei 
in the galaxy is 



to 



D 



10 7 yr 



R 



10 kpc 



26 10 19 eV 
~Z E 



2-<5 2 



(8) 



where R ~ 10 kpc is the diffusion distance scale. For dis- 
tances of the order of the Sun's Galactocentric distance, 
R ~ 8 kpc, this escape time is much longer than the time 
between individual GRBs, tD 3> tcRB ~ 10 4 — 10 5 yr, so 
the UHECR injection reaches a steady-state regime. 

We are interested in a steady-state solution, so that 
rii(E,r,t) should not depend on time. However, a few 
latest GRBs can cause fluctuations on the average, as 
discussed below. Let us first assume that Galactic halo 
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FIG. 2: Galactocentric anisotropy for a source distribution 
that traces the stellar counts in MW, modeled by random 
generation of 10 3 bursts separated by time intervals of 10 5 yr. 
The model parameters are the same as in Fig. 1. Although 
the anisotropy in protons is large at high energies, their con- 
tribution to the total flux is small, so the total anisotropy 
< 10%, consistent with the observations. The latest GRBs 
do not introduce a large degree of anisotropy, as it would be 
in the case of UHE protons, but they can create "hot spots" 
and clusters of events. 



is a sphere with radius Rq ~ 100 kpc and that all the 
sources are at the Galactic Center, so that the problem 
is spherically symmetric: Qi(E,f,t) — 5(r)Qo(E /E) 1 , 
rii(E,f,t) — rii(E,r). This is admittedly a simplified 
model, and we will replace it with a more realistic model 
below. Neglecting the energy losses inside the Galaxy, 
one obtains the solution of Eq. © with a boundary con- 
dition corresponding to a diminishing flux outside the 
Galaxy: 



rii(E,r) 



47rr Di{E) 



E 



(9) 



This solution corresponds to energy-dependent compo- 
sition for E > Eq. Indeed, at critical energy Eqj, which 
is different for each nucleus, the solution ([9]) changes 
from oc E~ 7 ~ Sl to oc E~i~ 2+S2 because of the change in 
Di(E). Since the change occurs at a rigidity-dependent 
critical energy Eo,i = eEoZi, the larger nuclei lag behind 
the lighter nuclei in terms of the critical energy and the 
change in slope. If protons dominate for E < Eq, their 
flux drops dramatically for E > Eq, and the heavier nu- 
clei dominate the flux. The higher Zi, the higher is the 
energy at which the species experiences a drop in flux. 

One can also understand the change in composition 
by considering the time of diffusion across the halo is 
ti ~ R 2 /Di. The longer the particle remains in the halo, 
the higher is the probability of its detection. At higher 
energies, the magnetic field's ability to delay the pas- 
sage of the particle diminishes, and the density of such 
particles drops precipitously for E > E o i . Since E^ is 
proportional to the electric charge, the drop in the flux 
occurs at different energies for different species. 



Of course, the assumption that all the sources are lo- 
cated in the Galactic center may not be realistic. If past 
GRBs in the Milky Way are the sources, one can model 
their distribution in different ways: one can assume (i) 
that all GRBs happen in the Galactic Center, or (ii) that 
GRB distribution follows the distribution of stars in MW, 
or (iii) one can include the short GRB distribution, which 
is expected to extend more into the halo. 

In Fig. 1 we show the spectrum calculated numerically 
for the source distribution (ii) , which we model using the 
star counts from Ref. [Hj]. Some 10 3 GRBs separated by 
time intervals of 10 5 years were generated in each Monte 
Carlo simulation, and the parameters were chosen to fit 
the data. We have assumed a two-component composi- 
tion with protons and iron nuclei. The best fit for 7 = 2.3 
is obtained for 90% protons and 10% iron, and » 4/iG 
magnetic field coherent on 0.2 kpc scale. For the case of 
short GRBs (iii), the distribution of sources can be ob- 
tained from observations [l9j]. The spectra obtained in 
this case are similar to those shown in Fig. 1. 

It is intriguing that the change from proton to iron in 
Fig. 1 is consistent with the dip in the spectrum that is 
usually attributed to either pair production or the change 
from Galactic to extragalactic component 0, [2l[ . How- 
ever, one should not consider the fit in Fig. 1 more than 
an illustration of the general principle. One must include 
multiple species of nuclei and the extragalactic protons, 
and one must model the propagation of UHECRs more 
carefully to compare the predictions with the data [l[ 
quantitatively. 

The data do not show a significant Galactocentric 
anisotropy in the arrival directions of UHECRs (although 
some clusters reported by PAO tend to gravitate toward 
the Galactic plane). For nuclei, however, one does not 
expect much anisotropy even if the sources are Galac- 
tic. We define the anisotropy parameter 8(E) in terms 
of maximal and minimal fluxes J m - m (E) and J mia c(E), 
depending on the arrival directions. In the diffusion ap- 
proximation, 



5(E) = 



"max J n 



d 



J: ti 



J n 



ZD(E)— In ni(E,r). (10) 



Obviously, model (i), assuming that all the sources are 
in the Galactic Center, predicts the largest anisotropy, 
hence setting the upper bound on 6. We find 5 < 0.1 for 
£<3x 10 19 eV. 

For model (ii), which assumes that the source dis- 
tribution follows the stellar distribution, the anisotropy 
can be computed numerically. We have calculated the 
anisotropy parameter by generating 10 3 GRBs occurring 
once every 10 5 years. The results are shown in Fig. 2 
for the same parameters as in Fig. 1. The anisotropy for 
model (iii) is even smaller. 

While the average flux includes contributions of GRBs 
form different distances and different times, the latest 
nearby GRBs can create fluctuations. A cluster of several 
UHECRs around Cen A detected by PAO may be the 
result of such a fluctuation due to one GRB that happens 
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to coincide with Cen A. Alternatively, since we expect the 
high-energy protons to escape from our Galaxy and from 
other galaxies, the cluster around Cen A may be due 
to extragalactic protons. Unlike protons, UHE nuclei 
from the last GRB do not introduce a large degree of 
anisotropy, as one can see form Fig. 2 based on a semi- 
realistic Monte Carlo simulation. 

Our model can be improved. First, one can use a more 
realistic source population model. Second, one should in- 
clude the coherent component of the Galactic magnetic 
field. Third, one should not assume that UHECRs com- 
prise only two types of particles, and one should include 
a realistic distribution of nuclei. Finally, one should in- 
clude the extragalactic component of UHECRs produced 
by distant sources, such as active galactic nuclei (AGN) 
and GRBs (outside MW) . A recent realization that very 
high energy gamma rays observed by Cherenkov tele- 
scopes from distant blazars are likely to be secondary 
photons produced in cosmic ray interactions along the 
line of sight lends further support to the assumption that 
cosmic rays are copiously produced in AGN jets [22| ■ For 
energies E > 3 x 10 19 eV, the energy losses due to pho- 
todisintegration, pion production, pair production and 
interactions with interstellar medium become important 



and must be included. The propagation distance in the 
Galaxy exceeds 10 Mpc, according to Eq. flSJ), so that 
the Galactic component should exhibit an analog of GZK 
suppression in the spectrum. The extragalactic propaga- 
tion can also affect the composition around 10 18 eV [2fJ, 
as well as the arrival direction anisotropy [231 ] . 

We have shown that effects of rigidity-dependent dif- 
fusion of UHECRs from possible Galactic sources, such 
as past GRBs (or rare types of supernovae), in the Milky 
Way can produce the energy-dependent composition as 
observed by PAO The simplest two-component 

model including protons and iron nuclei from the Galac- 
tic sources gives a good fit to composition [l| and spec- 
trum [24j for reasonable values of the Galactic magnetic 
fields. 
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